Modelamiento estadístico.
Primer congreso chileno de estudios interdisciplinarios sobre diversidad sexual y de género”

Autores/as
Afiliación

José Daniel Conejeros

Pontificia Universidad Católica de Chile

Javiera Porcel

Fecha de publicación

19 de marzo de 2025

¿Qué haremos en este bloque?

Cuando trabajamos con diseños complejos no podemos usar métodos estadísticos tradicionales sin ajustar la varianza y los ponderadores muestrales. En este bloque exploraremos modelos de regresión lineal, logística y técnicas avanzadas para análisis de encuestas.

Introducción a los modelos con encuestas

Preparamos los datos

Cargamos las librerías de trabajo:

Click para ver el código
# Para instalar manualmente 
install.packages("tidyverse")
install.packages("rio")
install.packages("haven")
install.packages("survey")
install.packages("srvyr")
install.packages("vtable")
install.packages("tidymodels")
install.packages("texreg")
install.packages("factoextra")
install.packages("NbClust")

# Cargamos las librerías de trabajo 
library("tidyverse")
library("rio")
library("haven")
library("survey")
library("srvyr")
library("vtable")
library("tidymodels")
library("texreg")
library("factoextra")
library("NbClust")

Abrimos los datos:

Click para ver el código
# Hay varios caminos para abrir los datos dependiendo el formato:
data <- rio::import("20240516_enssex_data.rdata") # R

Preparamos las variables:

Click para ver el código
variables <- c("folio_encuesta", "macro_region", "region", 
               "varstrat", # Pseuestrato
               "varunit", # Pseudoconglomerado
               "w_personas_cal", "w_personas_cal_nr_mod_hogar", # Ponderadores
               "fecha", "edad_primera_rsex", "relaciones_sexuales", 
               "p1", "genero_separada", "orientacion_sexual_separado",
               "edad", "fecha_nacimiento_ano", 
               "p34", # Conversaciones sobre sexualidad
               "i_6_p9", # Bienestar con la vida sexual
               "educ", 
               "p247", # Victima de abuso sexual
               "evaluacion_educacion_sexual", 
               "lugar_primera_rsex", 
               "vinculo_primera_rsex")

data <- data |> 
  dplyr::select(all_of(variables)) |> 
  dplyr::mutate(fecha_encuesta=dmy(fecha),
                agno_encuesta=year(fecha_encuesta), 
                edad_val=agno_encuesta-edad, 
                validacion_fecha_nac=if_else(edad_val==fecha_nacimiento_ano, 1, 0)) |> 
  mutate(
    macro_region=factor(macro_region, labels=c("Macrozona Norte",
                                               "Macrozona Centro", 
                                               "Macrozona Sur")),
    region=factor(region, labels=c("Tarapacá", 
                                   "Antofagasta", 
                                   "Atacama", 
                                   "Coquimbo",
                                   "Valparaíso", 
                                   "O'Higgins", 
                                   "Maule", 
                                   "Biobío", 
                                   "La Araucanía", 
                                   "Los Lagos", 
                                   "Aysén", 
                                   "Magallanes", 
                                   "Metropolitana", 
                                   "Los Ríos", 
                                   "Arica y Parinacota", 
                                   "Ñuble")),
    sex=factor(p1, labels=c("Hombre", "Mujer")), 
    genero=factor(genero_separada, labels=c("Masculino",
                                            "Transgénero", 
                                            "No binario", 
                                            "Otro",
                                            "Femenino",
                                            "No responder")),
    vinculo_primera_rsex=factor(vinculo_primera_rsex, labels = c("Pareja", 
                                                       "Amigo/a",
                                                       "Amigo/a con beneficios", 
                                                       "Otro",
                                                       "No responde")),
     lugar_primera_rsex=factor(lugar_primera_rsex, labels = c("Casa propia", 
                                                       "Casa de los padres",
                                                       "Casa de amigos/as", 
                                                       "Hotel",
                                                       "Aire libre",
                                                       "Otro",
                                                       "No responde")),
    cohorte=case_when(
      fecha_nacimiento_ano<1950 ~ 1, 
      fecha_nacimiento_ano>=1950 & fecha_nacimiento_ano<1960 ~ 2,
      fecha_nacimiento_ano>=1960 & fecha_nacimiento_ano<1970 ~ 3,
      fecha_nacimiento_ano>=1970 & fecha_nacimiento_ano<1980 ~ 4,
      fecha_nacimiento_ano>=1980 & fecha_nacimiento_ano<1990 ~ 5,
      fecha_nacimiento_ano>=1990 & fecha_nacimiento_ano<2000 ~ 6,
      fecha_nacimiento_ano>=2000 ~ 7,
      TRUE ~ NA
    ), 
    cohorte=factor(cohorte, labels = c("Cohorte < 1950", 
                                     "Cohorte 1951-1960",
                                     "Cohorte 1961-1970",
                                     "Cohorte 1971-1980",
                                     "Cohorte 1981-1990",
                                     "Cohorte 1991-2000",
                                     "Cohorte 2001-2005")),
    p34=if_else(p34 %in% c(8,9), 3, p34-1),
    p34=if_else(p34==2, 1, p34),
    sexual_convers=factor(p34, labels = c("No", "Si", "No responde")), 
    educ = factor(educ, labels=c("Primaria",
                                 "Secundaria incompleta",
                                 "Secondarian completa",
                                 "Universitario")),
    sexual_abuso = factor(p247, labels=c("Si", "No", "No responde")),
    sexual_educ=factor(evaluacion_educacion_sexual, labels=c("Mal", 
                                                             "Regular", 
                                                             "Bien", 
                                                             "Sin respuesta"
                                                                  )),

    bienestar=if_else(i_6_p9==9, NA, i_6_p9),
    bienestar=as.numeric(bienestar)
  ) |> 
  dplyr::select(varstrat, varunit, w_personas_cal, 
         macro_region, region,
         bienestar,
         edad_primera_rsex, 
         vinculo_primera_rsex, lugar_primera_rsex,
         sex, orientacion_sexual_separado, 
         edad, cohorte, educ,
         sexual_convers,
         sexual_educ, 
         sexual_abuso
         ) |> 
  drop_na() 

Revisamos nuestra tabla descriptiva:

Click para ver el código
# Tabla de descriptivos sin ponderar 
st(dplyr::select(data, !c("varstrat", "varunit", "w_personas_cal", "region")),
   #labels=c("Identidad de género", "Grupo Etario", "Nivel educativo", "Víctima de violencia"),
   digits = 3, 
   out="kable", 
   add.median = TRUE,
   fixed.digits = TRUE, 
   simple.kable = FALSE,
   title="",
   numformat = NA) %>%
  kableExtra::kable_styling(latex_options = c("scale_down", "hold_position"),
                            full_width = FALSE, fixed_thead = T)
Variable N Mean Std. Dev. Min Pctl. 25 Pctl. 50 Pctl. 75 Max
macro_region 16904
... Macrozona Norte 3261 19.3%
... Macrozona Centro 8457 50.0%
... Macrozona Sur 5186 30.7%
bienestar 16904 4.964 1.997 1 4 6 7 7
edad_primera_rsex 16904 17.739 3.484 10 16 17 19 76
vinculo_primera_rsex 16904
... Pareja 10319 61.0%
... Amigo/a 2164 12.8%
... Amigo/a con beneficios 2267 13.4%
... Otro 2082 12.3%
... No responde 72 0.4%
lugar_primera_rsex 16904
... Casa propia 4812 28.5%
... Casa de los padres 5434 32.1%
... Casa de amigos/as 1686 10.0%
... Hotel 893 5.3%
... Aire libre 1893 11.2%
... Otro 1824 10.8%
... No responde 362 2.1%
sex 16904
... Hombre 5805 34.3%
... Mujer 11099 65.7%
edad 16904 45.104 17.123 18 30 43 58 100
cohorte 16904
... Cohorte < 1950 1172 6.9%
... Cohorte 1951-1960 1986 11.7%
... Cohorte 1961-1970 2693 15.9%
... Cohorte 1971-1980 2806 16.6%
... Cohorte 1981-1990 3139 18.6%
... Cohorte 1991-2000 3814 22.6%
... Cohorte 2001-2005 1294 7.7%
educ 16904
... Primaria 3033 17.9%
... Secundaria incompleta 2365 14.0%
... Secondarian completa 5765 34.1%
... Universitario 5741 34.0%
sexual_convers 16904
... No 11893 70.4%
... Si 4841 28.6%
... No responde 170 1.0%
sexual_educ 16904
... Mal 8167 48.3%
... Regular 4311 25.5%
... Bien 3571 21.1%
... Sin respuesta 855 5.1%
sexual_abuso 16904
... Si 2106 12.5%
... No 13647 80.7%
... No responde 1151 6.8%

Regresión lineal

La regresión lineal ponderada estima la relación entre una variable dependiente continua \(Y\) y un conjunto de predictores \(X_1, X_2, \dots, X_p\):

\[Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \dots + \beta_p X_{ip} + \varepsilon_i\]

Donde:
- \(Y_i\) es la variable dependiente (Ej: edad de inicio sexual).
- \(X_{i1}, X_{i2}, \dots, X_{ip}\) son las variables explicativas (Ej: sexo, nivel educativo, región).
- \(\beta_0\) es la intersección (constante).
- \(\beta_p\) son los coeficientes de regresión.
- \(\varepsilon_i\) es el error aleatorio.

Click para ver el código
# Sin diseño muestral
m1 <- lm(bienestar ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data=data)

# Podemos ver los resultados de distintas maneras:
summary(m1)

Call:
lm(formula = bienestar ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
    edad + cohorte + educ + sexual_convers + sexual_educ, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-5.231 -1.158  0.423  1.446  4.492 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  5.684779   0.396137  14.351  < 2e-16 ***
edad_primera_rsex           -0.017465   0.004425  -3.947 7.96e-05 ***
sexMujer                    -0.360791   0.030858 -11.692  < 2e-16 ***
orientacion_sexual_separado -0.052982   0.012477  -4.246 2.18e-05 ***
edad                        -0.026080   0.004952  -5.267 1.40e-07 ***
cohorteCohorte 1951-1960     0.502219   0.087169   5.761 8.49e-09 ***
cohorteCohorte 1961-1970     0.731319   0.120681   6.060 1.39e-09 ***
cohorteCohorte 1971-1980     0.861353   0.162963   5.286 1.27e-07 ***
cohorteCohorte 1981-1990     0.831484   0.210517   3.950 7.86e-05 ***
cohorteCohorte 1991-2000     0.597319   0.255378   2.339 0.019349 *  
cohorteCohorte 2001-2005     0.271129   0.294317   0.921 0.356951    
educSecundaria incompleta    0.170301   0.051936   3.279 0.001044 ** 
educSecondarian completa     0.388808   0.045655   8.516  < 2e-16 ***
educUniversitario            0.464984   0.047432   9.803  < 2e-16 ***
sexual_conversSi             0.091151   0.034890   2.613 0.008995 ** 
sexual_conversNo responde    0.032265   0.144041   0.224 0.822760    
sexual_educRegular           0.118845   0.035977   3.303 0.000957 ***
sexual_educBien              0.296057   0.039100   7.572 3.87e-14 ***
sexual_educSin respuesta     0.117389   0.067402   1.742 0.081592 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.856 on 16885 degrees of freedom
Multiple R-squared:  0.1374,    Adjusted R-squared:  0.1365 
F-statistic: 149.4 on 18 and 16885 DF,  p-value: < 2.2e-16
Click para ver el código
tidy(m1)
# A tibble: 19 × 5
   term                        estimate std.error statistic  p.value
   <chr>                          <dbl>     <dbl>     <dbl>    <dbl>
 1 (Intercept)                   5.68     0.396      14.4   1.98e-46
 2 edad_primera_rsex            -0.0175   0.00443    -3.95  7.96e- 5
 3 sexMujer                     -0.361    0.0309    -11.7   1.86e-31
 4 orientacion_sexual_separado  -0.0530   0.0125     -4.25  2.18e- 5
 5 edad                         -0.0261   0.00495    -5.27  1.40e- 7
 6 cohorteCohorte 1951-1960      0.502    0.0872      5.76  8.49e- 9
 7 cohorteCohorte 1961-1970      0.731    0.121       6.06  1.39e- 9
 8 cohorteCohorte 1971-1980      0.861    0.163       5.29  1.27e- 7
 9 cohorteCohorte 1981-1990      0.831    0.211       3.95  7.86e- 5
10 cohorteCohorte 1991-2000      0.597    0.255       2.34  1.93e- 2
11 cohorteCohorte 2001-2005      0.271    0.294       0.921 3.57e- 1
12 educSecundaria incompleta     0.170    0.0519      3.28  1.04e- 3
13 educSecondarian completa      0.389    0.0457      8.52  1.78e-17
14 educUniversitario             0.465    0.0474      9.80  1.25e-22
15 sexual_conversSi              0.0912   0.0349      2.61  8.99e- 3
16 sexual_conversNo responde     0.0323   0.144       0.224 8.23e- 1
17 sexual_educRegular            0.119    0.0360      3.30  9.57e- 4
18 sexual_educBien               0.296    0.0391      7.57  3.87e-14
19 sexual_educSin respuesta      0.117    0.0674      1.74  8.16e- 2
Click para ver el código
screenreg(l=list(m1), digits = 3)

==========================================
                             Model 1      
------------------------------------------
(Intercept)                      5.685 ***
                                (0.396)   
edad_primera_rsex               -0.017 ***
                                (0.004)   
sexMujer                        -0.361 ***
                                (0.031)   
orientacion_sexual_separado     -0.053 ***
                                (0.012)   
edad                            -0.026 ***
                                (0.005)   
cohorteCohorte 1951-1960         0.502 ***
                                (0.087)   
cohorteCohorte 1961-1970         0.731 ***
                                (0.121)   
cohorteCohorte 1971-1980         0.861 ***
                                (0.163)   
cohorteCohorte 1981-1990         0.831 ***
                                (0.211)   
cohorteCohorte 1991-2000         0.597 *  
                                (0.255)   
cohorteCohorte 2001-2005         0.271    
                                (0.294)   
educSecundaria incompleta        0.170 ** 
                                (0.052)   
educSecondarian completa         0.389 ***
                                (0.046)   
educUniversitario                0.465 ***
                                (0.047)   
sexual_conversSi                 0.091 ** 
                                (0.035)   
sexual_conversNo responde        0.032    
                                (0.144)   
sexual_educRegular               0.119 ***
                                (0.036)   
sexual_educBien                  0.296 ***
                                (0.039)   
sexual_educSin respuesta         0.117    
                                (0.067)   
------------------------------------------
R^2                              0.137    
Adj. R^2                         0.137    
Num. obs.                    16904        
==========================================
*** p < 0.001; ** p < 0.01; * p < 0.05
Click para ver el código
# Con diseño muestral
# Definir diseño muestral ENSSEX
enssex_design <- svydesign(
  id = ~varunit,   # UPM
  strata = ~varstrat,    # Estratos
  weights = ~ w_personas_cal,  # Ponderador muestral
  data = data,
  nest = TRUE,
)

options(survey.lonely.psu = "adjust")

# Estimamos
m2 <- svyglm(bienestar ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, design = enssex_design)

screenreg(l=list(m1, m2), digits = 3)

=========================================================
                             Model 1        Model 2      
---------------------------------------------------------
(Intercept)                      5.685 ***      5.387 ***
                                (0.396)        (0.554)   
edad_primera_rsex               -0.017 ***     -0.030 ***
                                (0.004)        (0.007)   
sexMujer                        -0.361 ***     -0.409 ***
                                (0.031)        (0.048)   
orientacion_sexual_separado     -0.053 ***     -0.063 ** 
                                (0.012)        (0.019)   
edad                            -0.026 ***     -0.017 *  
                                (0.005)        (0.007)   
cohorteCohorte 1951-1960         0.502 ***      0.554 ***
                                (0.087)        (0.115)   
cohorteCohorte 1961-1970         0.731 ***      0.911 ***
                                (0.121)        (0.165)   
cohorteCohorte 1971-1980         0.861 ***      1.036 ***
                                (0.163)        (0.225)   
cohorteCohorte 1981-1990         0.831 ***      1.015 ***
                                (0.211)        (0.281)   
cohorteCohorte 1991-2000         0.597 *        0.814 *  
                                (0.255)        (0.344)   
cohorteCohorte 2001-2005         0.271          0.631    
                                (0.294)        (0.400)   
educSecundaria incompleta        0.170 **       0.268 ** 
                                (0.052)        (0.088)   
educSecondarian completa         0.389 ***      0.526 ***
                                (0.046)        (0.076)   
educUniversitario                0.465 ***      0.579 ***
                                (0.047)        (0.079)   
sexual_conversSi                 0.091 **       0.087 *  
                                (0.035)        (0.043)   
sexual_conversNo responde        0.032         -0.330    
                                (0.144)        (0.191)   
sexual_educRegular               0.119 ***      0.129 ** 
                                (0.036)        (0.049)   
sexual_educBien                  0.296 ***      0.279 ***
                                (0.039)        (0.058)   
sexual_educSin respuesta         0.117         -0.057    
                                (0.067)        (0.094)   
---------------------------------------------------------
R^2                              0.137                   
Adj. R^2                         0.137                   
Num. obs.                    16904          16904        
Deviance                                    51158.348    
Dispersion                                      3.027    
=========================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Regresión logística

Cuando la variable dependiente es binaria \(Y \in {0,1}\), usamos regresión logística ponderada:

Donde: - \(P(Y_i = 1)\) es la probabilidad de que ocurra el evento (Ej: uso de anticonceptivos = 1, no uso = 0). - \(\frac{P(Y_i = 1)}{1 - P(Y_i = 1)}\) es el odds ratio (razón de probabilidades). - \(\beta_p\) representa el efecto de cada variable independiente sobre la probabilidad del evento.

Click para ver el código
# Modelo sin diseño muestral
l1 <- glm(factor(bienestar>6) ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ,
        data=data,
        family = binomial(link="logit"))

# Modelo con diseño muestral
l2 <- svyglm(factor(bienestar>6) ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ,
       design = enssex_design, 
       family = quasibinomial())

# Reportamos
screenreg(l=list(l1, l2), digits = 3)

=========================================================
                             Model 1        Model 2      
---------------------------------------------------------
(Intercept)                     -0.595         -1.543 *  
                                (0.506)        (0.739)   
edad_primera_rsex               -0.019 **      -0.028 ** 
                                (0.006)        (0.009)   
sexMujer                        -0.017         -0.059    
                                (0.038)        (0.065)   
orientacion_sexual_separado     -0.030         -0.029    
                                (0.016)        (0.030)   
edad                            -0.016 *       -0.003    
                                (0.006)        (0.009)   
cohorteCohorte 1951-1960         0.286 *        0.448 *  
                                (0.133)        (0.197)   
cohorteCohorte 1961-1970         0.547 **       0.940 ***
                                (0.167)        (0.232)   
cohorteCohorte 1971-1980         0.711 ***      1.212 ***
                                (0.216)        (0.307)   
cohorteCohorte 1981-1990         0.720 **       1.229 ** 
                                (0.274)        (0.385)   
cohorteCohorte 1991-2000         0.641          1.236 ** 
                                (0.329)        (0.464)   
cohorteCohorte 2001-2005         0.509          1.248 *  
                                (0.377)        (0.535)   
educSecundaria incompleta        0.006          0.088    
                                (0.069)        (0.115)   
educSecondarian completa         0.063          0.213 *  
                                (0.059)        (0.099)   
educUniversitario                0.085          0.214 *  
                                (0.061)        (0.102)   
sexual_conversSi                 0.111 **       0.123    
                                (0.041)        (0.064)   
sexual_conversNo responde       -0.197         -0.534    
                                (0.179)        (0.296)   
sexual_educRegular              -0.031          0.029    
                                (0.044)        (0.065)   
sexual_educBien                  0.299 ***      0.384 ***
                                (0.046)        (0.079)   
sexual_educSin respuesta         0.160         -0.077    
                                (0.088)        (0.122)   
---------------------------------------------------------
AIC                          19346.677                   
BIC                          19493.647                   
Log Likelihood               -9654.338                   
Deviance                     19308.677      19842.472    
Num. obs.                    16904          16904        
Dispersion                                      1.001    
=========================================================
*** p < 0.001; ** p < 0.01; * p < 0.05
Click para ver el código
tidy(l1, exponentiate = TRUE)
# A tibble: 19 × 5
   term                        estimate std.error statistic  p.value
   <chr>                          <dbl>     <dbl>     <dbl>    <dbl>
 1 (Intercept)                    0.552   0.506     -1.17   2.40e- 1
 2 edad_primera_rsex              0.981   0.00597   -3.27   1.09e- 3
 3 sexMujer                       0.983   0.0376    -0.445  6.56e- 1
 4 orientacion_sexual_separado    0.970   0.0157    -1.93   5.32e- 2
 5 edad                           0.984   0.00628   -2.52   1.18e- 2
 6 cohorteCohorte 1951-1960       1.33    0.133      2.15   3.19e- 2
 7 cohorteCohorte 1961-1970       1.73    0.167      3.28   1.04e- 3
 8 cohorteCohorte 1971-1980       2.04    0.216      3.29   9.91e- 4
 9 cohorteCohorte 1981-1990       2.05    0.274      2.63   8.56e- 3
10 cohorteCohorte 1991-2000       1.90    0.329      1.95   5.14e- 2
11 cohorteCohorte 2001-2005       1.66    0.377      1.35   1.76e- 1
12 educSecundaria incompleta      1.01    0.0688     0.0827 9.34e- 1
13 educSecondarian completa       1.06    0.0594     1.06   2.90e- 1
14 educUniversitario              1.09    0.0610     1.39   1.63e- 1
15 sexual_conversSi               1.12    0.0409     2.71   6.67e- 3
16 sexual_conversNo responde      0.821   0.179     -1.10   2.71e- 1
17 sexual_educRegular             0.970   0.0441    -0.700  4.84e- 1
18 sexual_educBien                1.35    0.0460     6.51   7.67e-11
19 sexual_educSin respuesta       1.17    0.0881     1.82   6.91e- 2
Click para ver el código
tidy(l2, exponentiate = TRUE)
# A tibble: 19 × 5
   term                        estimate std.error statistic    p.value
   <chr>                          <dbl>     <dbl>     <dbl>      <dbl>
 1 (Intercept)                    0.214   0.739      -2.09  0.0377    
 2 edad_primera_rsex              0.973   0.00935    -2.94  0.00351   
 3 sexMujer                       0.942   0.0646     -0.917 0.360     
 4 orientacion_sexual_separado    0.971   0.0297     -0.981 0.328     
 5 edad                           0.997   0.00894    -0.315 0.753     
 6 cohorteCohorte 1951-1960       1.57    0.197       2.27  0.0238    
 7 cohorteCohorte 1961-1970       2.56    0.232       4.05  0.0000671 
 8 cohorteCohorte 1971-1980       3.36    0.307       3.95  0.0000990 
 9 cohorteCohorte 1981-1990       3.42    0.385       3.19  0.00158   
10 cohorteCohorte 1991-2000       3.44    0.464       2.67  0.00813   
11 cohorteCohorte 2001-2005       3.48    0.535       2.33  0.0205    
12 educSecundaria incompleta      1.09    0.115       0.765 0.445     
13 educSecondarian completa       1.24    0.0985      2.16  0.0313    
14 educUniversitario              1.24    0.102       2.09  0.0374    
15 sexual_conversSi               1.13    0.0637      1.93  0.0543    
16 sexual_conversNo responde      0.586   0.296      -1.80  0.0726    
17 sexual_educRegular             1.03    0.0653      0.450 0.653     
18 sexual_educBien                1.47    0.0788      4.87  0.00000188
19 sexual_educSin respuesta       0.926   0.122      -0.634 0.527     

Aplicación I: análisis de sobrevivencia

El modelo de Cox (o regresión de riesgos proporcionales) se usa para tiempo hasta un evento (Ej: edad de primera relación sexual). Se expresa como:

\[h(t | X) = h_0(t) \exp(\beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p)\]

Donde: - \(h(t | X)\) es el riesgo instantáneo de que ocurra el evento en el tiempo ( t ). - \(h_0(t)\) es la función de riesgo base. - \(\beta_p\) son los coeficientes de las variables explicativas ( X_p ).

Click para ver el código
# Modelos de COX
cox <- coxph(Surv(edad_primera_rsex) ~ bienestar + sex + orientacion_sexual_separado + edad + cohorte + educ + sexual_convers + sexual_educ, data = data)

cox_ponderado <- svycoxph(Surv(edad_primera_rsex) ~ bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, design = enssex_design)

# Modelos de AFT
aft1 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "gaussian")
aft2 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "logistic")
aft3 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "extreme")
aft4 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "weibull")
aft5 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "exponential")
aft6 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "lognormal")

# Reportamos 
summary(cox)
Call:
coxph(formula = Surv(edad_primera_rsex) ~ bienestar + sex + orientacion_sexual_separado + 
    edad + cohorte + educ + sexual_convers + sexual_educ, data = data)

  n= 16904, number of events= 16904 

                                 coef exp(coef)  se(coef)       z Pr(>|z|)    
bienestar                    0.017808  1.017967  0.004111   4.332 1.48e-05 ***
sexMujer                    -0.388540  0.678046  0.016454 -23.614  < 2e-16 ***
orientacion_sexual_separado  0.011343  1.011408  0.006704   1.692  0.09065 .  
edad                        -0.024532  0.975766  0.002706  -9.065  < 2e-16 ***
cohorteCohorte 1951-1960    -0.134674  0.874001  0.046993  -2.866  0.00416 ** 
cohorteCohorte 1961-1970    -0.314725  0.729990  0.065426  -4.810 1.51e-06 ***
cohorteCohorte 1971-1980    -0.260730  0.770489  0.088389  -2.950  0.00318 ** 
cohorteCohorte 1981-1990    -0.260043  0.771018  0.114083  -2.279  0.02264 *  
cohorteCohorte 1991-2000    -0.252746  0.776665  0.138491  -1.825  0.06800 .  
cohorteCohorte 2001-2005    -0.107290  0.898265  0.159751  -0.672  0.50183    
educSecundaria incompleta   -0.077818  0.925133  0.028045  -2.775  0.00552 ** 
educSecondarian completa    -0.180989  0.834445  0.024533  -7.377 1.61e-13 ***
educUniversitario           -0.314788  0.729944  0.025438 -12.375  < 2e-16 ***
sexual_conversSi            -0.015008  0.985104  0.018770  -0.800  0.42395    
sexual_conversNo responde    0.007278  1.007305  0.077827   0.094  0.92549    
sexual_educRegular          -0.038835  0.961910  0.019435  -1.998  0.04570 *  
sexual_educBien             -0.059724  0.942024  0.021169  -2.821  0.00478 ** 
sexual_educSin respuesta    -0.028206  0.972188  0.036386  -0.775  0.43823    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                            exp(coef) exp(-coef) lower .95 upper .95
bienestar                      1.0180     0.9823    1.0098    1.0262
sexMujer                       0.6780     1.4748    0.6565    0.7003
orientacion_sexual_separado    1.0114     0.9887    0.9982    1.0248
edad                           0.9758     1.0248    0.9706    0.9810
cohorteCohorte 1951-1960       0.8740     1.1442    0.7971    0.9583
cohorteCohorte 1961-1970       0.7300     1.3699    0.6421    0.8299
cohorteCohorte 1971-1980       0.7705     1.2979    0.6479    0.9162
cohorteCohorte 1981-1990       0.7710     1.2970    0.6165    0.9642
cohorteCohorte 1991-2000       0.7767     1.2876    0.5920    1.0189
cohorteCohorte 2001-2005       0.8983     1.1133    0.6568    1.2285
educSecundaria incompleta      0.9251     1.0809    0.8757    0.9774
educSecondarian completa       0.8344     1.1984    0.7953    0.8755
educUniversitario              0.7299     1.3700    0.6944    0.7673
sexual_conversSi               0.9851     1.0151    0.9495    1.0220
sexual_conversNo responde      1.0073     0.9927    0.8648    1.1733
sexual_educRegular             0.9619     1.0396    0.9260    0.9993
sexual_educBien                0.9420     1.0615    0.9037    0.9819
sexual_educSin respuesta       0.9722     1.0286    0.9053    1.0441

Concordance= 0.647  (se = 0.003 )
Likelihood ratio test= 2755  on 18 df,   p=<2e-16
Wald test            = 2838  on 18 df,   p=<2e-16
Score (logrank) test = 2939  on 18 df,   p=<2e-16
Click para ver el código
summary(cox_ponderado)
Stratified 1 - level Cluster Sampling design (with replacement)
With (372) clusters.
svydesign(id = ~varunit, strata = ~varstrat, weights = ~w_personas_cal, 
    data = data, nest = TRUE, )
Call:
svycoxph(formula = Surv(edad_primera_rsex) ~ bienestar + sex + 
    orientacion_sexual_separado + edad + cohorte + educ + sexual_convers + 
    sexual_educ, design = enssex_design)

  n= 16904, number of events= 16904 

                                  coef  exp(coef)   se(coef)  robust se       z
bienestar                    3.159e-02  1.032e+00  4.331e-03  6.418e-03   4.923
sexMujer                    -3.512e-01  7.039e-01  1.571e-02  2.544e-02 -13.803
orientacion_sexual_separado -4.932e-03  9.951e-01  8.071e-03  1.824e-02  -0.270
edad                        -2.475e-02  9.756e-01  2.696e-03  4.252e-03  -5.819
cohorteCohorte 1951-1960    -1.444e-01  8.655e-01  4.943e-02  7.105e-02  -2.033
cohorteCohorte 1961-1970    -3.588e-01  6.985e-01  6.747e-02  9.497e-02  -3.778
cohorteCohorte 1971-1980    -2.749e-01  7.597e-01  9.050e-02  1.379e-01  -1.994
cohorteCohorte 1981-1990    -3.729e-01  6.887e-01  1.152e-01  1.830e-01  -2.038
cohorteCohorte 1991-2000    -4.382e-01  6.452e-01  1.397e-01  2.146e-01  -2.041
cohorteCohorte 2001-2005    -2.968e-01  7.432e-01  1.612e-01  2.522e-01  -1.177
educSecundaria incompleta   -9.524e-02  9.092e-01  3.231e-02  4.782e-02  -1.992
educSecondarian completa    -2.573e-01  7.731e-01  2.617e-02  5.084e-02  -5.062
educUniversitario           -3.366e-01  7.142e-01  2.614e-02  4.758e-02  -7.074
sexual_conversSi            -5.884e-05  9.999e-01  1.852e-02  2.960e-02  -0.002
sexual_conversNo responde   -2.425e-02  9.760e-01  8.717e-02  1.090e-01  -0.223
sexual_educRegular          -2.107e-02  9.791e-01  1.951e-02  3.492e-02  -0.603
sexual_educBien              1.586e-02  1.016e+00  2.102e-02  3.351e-02   0.473
sexual_educSin respuesta    -7.453e-02  9.282e-01  3.196e-02  5.383e-02  -1.384
                            Pr(>|z|)    
bienestar                   8.52e-07 ***
sexMujer                     < 2e-16 ***
orientacion_sexual_separado 0.786901    
edad                        5.91e-09 ***
cohorteCohorte 1951-1960    0.042092 *  
cohorteCohorte 1961-1970    0.000158 ***
cohorteCohorte 1971-1980    0.046185 *  
cohorteCohorte 1981-1990    0.041569 *  
cohorteCohorte 1991-2000    0.041202 *  
cohorteCohorte 2001-2005    0.239312    
educSecundaria incompleta   0.046404 *  
educSecondarian completa    4.15e-07 ***
educUniversitario           1.51e-12 ***
sexual_conversSi            0.998414    
sexual_conversNo responde   0.823862    
sexual_educRegular          0.546203    
sexual_educBien             0.636021    
sexual_educSin respuesta    0.166217    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                            exp(coef) exp(-coef) lower .95 upper .95
bienestar                      1.0321     0.9689    1.0192    1.0452
sexMujer                       0.7039     1.4207    0.6696    0.7399
orientacion_sexual_separado    0.9951     1.0049    0.9601    1.0313
edad                           0.9756     1.0251    0.9675    0.9837
cohorteCohorte 1951-1960       0.8655     1.1554    0.7530    0.9949
cohorteCohorte 1961-1970       0.6985     1.4316    0.5799    0.8414
cohorteCohorte 1971-1980       0.7597     1.3164    0.5798    0.9954
cohorteCohorte 1981-1990       0.6887     1.4520    0.4811    0.9859
cohorteCohorte 1991-2000       0.6452     1.5499    0.4236    0.9827
cohorteCohorte 2001-2005       0.7432     1.3455    0.4533    1.2184
educSecundaria incompleta      0.9092     1.0999    0.8278    0.9985
educSecondarian completa       0.7731     1.2935    0.6998    0.8541
educUniversitario              0.7142     1.4001    0.6506    0.7840
sexual_conversSi               0.9999     1.0001    0.9436    1.0597
sexual_conversNo responde      0.9760     1.0245    0.7884    1.2084
sexual_educRegular             0.9791     1.0213    0.9144    1.0485
sexual_educBien                1.0160     0.9843    0.9514    1.0850
sexual_educSin respuesta       0.9282     1.0774    0.8352    1.0315

Concordance= 0.629  (se = 0.004 )
Likelihood ratio test= NA  on 18 df,   p=NA
Wald test            = 1197  on 18 df,   p=<2e-16
Score (logrank) test = NA  on 18 df,   p=NA

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).
Click para ver el código
screenreg(l=list(cox, aft1, aft2, aft3, aft4, aft5, aft6), digits = 3,
        custom.model.names = c("Cox", "AFT Gaussiano", "AFT Logístico", "AFT Extremo", 
                               "AFT Weibull", "AFT Exponential", "AFT lognormal"))

============================================================================================================================================
                             Cox             AFT Gaussiano   AFT Logístico   AFT Extremo     AFT Weibull     AFT Exponential  AFT lognormal 
--------------------------------------------------------------------------------------------------------------------------------------------
bienestar                         0.018 ***      -0.053 ***      -0.050 ***      -0.082 ***      -0.006 ***      -0.003           -0.003 ***
                                 (0.004)         (0.013)         (0.011)         (0.022)         (0.001)         (0.004)          (0.001)   
sexMujer                         -0.389 ***       1.449 ***       1.264 ***       2.076 ***       0.087 ***       0.080 ***        0.080 ***
                                 (0.016)         (0.053)         (0.044)         (0.088)         (0.003)         (0.016)          (0.003)   
orientacion_sexual_separado       0.011          -0.045 *        -0.033          -0.046          -0.003 *        -0.002           -0.002 *  
                                 (0.007)         (0.022)         (0.018)         (0.037)         (0.001)         (0.007)          (0.001)   
edad                             -0.025 ***       0.077 ***       0.063 ***       0.250 ***       0.007 ***       0.004            0.004 ***
                                 (0.003)         (0.009)         (0.007)         (0.013)         (0.001)         (0.003)          (0.000)   
cohorteCohorte 1951-1960         -0.135 **        0.305 *         0.258           2.100 ***       0.008           0.018            0.019 *  
                                 (0.047)         (0.152)         (0.133)         (0.258)         (0.010)         (0.047)          (0.008)   
cohorteCohorte 1961-1970         -0.315 ***       0.738 ***       0.630 ***       3.835 ***       0.042 **        0.043            0.042 ***
                                 (0.065)         (0.210)         (0.181)         (0.346)         (0.014)         (0.065)          (0.011)   
cohorteCohorte 1971-1980         -0.261 **        0.432           0.284           4.483 ***       0.042 *         0.028            0.026    
                                 (0.088)         (0.283)         (0.242)         (0.460)         (0.019)         (0.088)          (0.014)   
cohorteCohorte 1981-1990         -0.260 *         0.451           0.301           5.053 ***       0.042           0.028            0.026    
                                 (0.114)         (0.366)         (0.312)         (0.591)         (0.024)         (0.114)          (0.018)   
cohorteCohorte 1991-2000         -0.253           0.574           0.366           6.534 ***       0.060 *         0.033            0.029    
                                 (0.138)         (0.444)         (0.378)         (0.712)         (0.029)         (0.138)          (0.022)   
cohorteCohorte 2001-2005         -0.107           0.655           0.449           7.526 ***       0.061           0.034            0.031    
                                 (0.160)         (0.511)         (0.434)         (0.816)         (0.034)         (0.159)          (0.026)   
educSecundaria incompleta        -0.078 **        0.437 ***       0.477 ***       0.364 *         0.010           0.023            0.024 ***
                                 (0.028)         (0.090)         (0.078)         (0.151)         (0.006)         (0.028)          (0.005)   
educSecondarian completa         -0.181 ***       0.844 ***       0.864 ***       0.873 ***       0.023 ***       0.045            0.047 ***
                                 (0.025)         (0.079)         (0.069)         (0.131)         (0.005)         (0.025)          (0.004)   
educUniversitario                -0.315 ***       1.259 ***       1.271 ***       0.949 ***       0.043 ***       0.068 **         0.071 ***
                                 (0.025)         (0.082)         (0.071)         (0.136)         (0.005)         (0.025)          (0.004)   
sexual_conversSi                 -0.015           0.067           0.012           0.576 ***       0.007           0.003            0.003    
                                 (0.019)         (0.061)         (0.050)         (0.099)         (0.004)         (0.019)          (0.003)   
sexual_conversNo responde         0.007          -0.000           0.095          -0.467          -0.018           0.001            0.003    
                                 (0.078)         (0.250)         (0.205)         (0.412)         (0.016)         (0.078)          (0.013)   
sexual_educRegular               -0.039 *         0.162 **        0.146 **        0.609 ***       0.020 ***       0.009            0.008 ** 
                                 (0.019)         (0.063)         (0.053)         (0.103)         (0.004)         (0.019)          (0.003)   
sexual_educBien                  -0.060 **        0.182 **        0.187 **        0.274 *         0.031 ***       0.010            0.010 ** 
                                 (0.021)         (0.068)         (0.057)         (0.113)         (0.004)         (0.021)          (0.003)   
sexual_educSin respuesta         -0.028           0.115          -0.031           1.514 ***       0.056 ***       0.006            0.002    
                                 (0.036)         (0.117)         (0.102)         (0.196)         (0.008)         (0.036)          (0.006)   
(Intercept)                                      12.279 ***      12.863 ***       1.677           2.540 ***       2.573 ***        2.574 ***
                                                 (0.686)         (0.584)         (1.074)         (0.045)         (0.213)          (0.035)   
Log(scale)                                        1.171 ***       0.465 ***       1.666 ***      -1.560 ***                       -1.816 ***
                                                 (0.005)         (0.007)         (0.004)         (0.005)                          (0.005)   
--------------------------------------------------------------------------------------------------------------------------------------------
AIC                          292616.244       87594.104       84196.554      102673.920       91600.274      130979.129        83268.230    
R^2                               0.150                                                                                                     
Max. R^2                          1.000                                                                                                     
Num. events                   16904                                                                                                         
Num. obs.                     16904           16904           16904           16904           16904           16904            16904        
Missings                          0                                                                                                         
PH test                           0.000                                                                                                     
BIC                                           87748.810       84351.260      102828.626       91754.980      131126.100        83422.936    
Log Likelihood                               -43777.052      -42078.277      -51316.960      -45780.137      -65470.564       -41614.115    
============================================================================================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Aplicación II: Análisis de clusters y segmentación

Agrupa observaciones según similitudes en sus características:

\[\arg\min_S \sum_{i=1}^{k} \sum_{x \in S_i} | x - \mu_i |^2\]

Donde: - \(S_i\) son los grupos (clusters). - \(\mu_i\) es el centroide de cada grupo. - \(x\) son los valores de las observaciones.

Para esto utilizaremos un set de preguntas:

Click para ver el código
# Usemos un set de preguntas
data <- rio::import("https://datos.gob.cl/dataset/c6983439-49f6-4e71-85fe-e8de6e73dae0/resource/ed81f50c-1c7d-43d9-9083-dfc161e0cd66/download/20240516_enssex_data.rdata") |> 
  select(ends_with("p9")) |> 
  mutate(across(everything(), ~ zap_labels(na_if(., 9)))) |> 
  drop_na()
Click para ver el código
glimpse(data)
Rows: 18,995
Columns: 6
$ i_1_p9 <dbl> 7, 4, 4, 6, 1, 4, 3, 7, 7, 1, 1, 3, 1, 3, 7, 7, 4, 7, 5, 5, 4, …
$ i_2_p9 <dbl> 6, 4, 4, 4, 1, 5, 7, 7, 7, 4, 4, 4, 3, 3, 5, 7, 6, 7, 6, 7, 3, …
$ i_3_p9 <dbl> 7, 3, 5, 5, 1, 6, 5, 7, 5, 2, 7, 2, 4, 4, 6, 7, 6, 7, 4, 5, 6, …
$ i_4_p9 <dbl> 6, 6, 4, 6, 3, 6, 5, 6, 7, 6, 1, 1, 7, 4, 5, 5, 3, 7, 5, 5, 6, …
$ i_5_p9 <dbl> 7, 6, 5, 4, 2, 5, 5, 6, 7, 7, 3, 6, 6, 5, 6, 7, 5, 6, 5, 5, 7, …
$ i_6_p9 <dbl> 7, 6, 5, 4, 6, 4, 5, 6, 5, 1, 4, 2, 6, 5, 6, 7, 6, 7, 4, 5, 4, …

Veamos los descriptivos (sin ponderar):

Click para ver el código
st(data,
   digits = 3, 
   out="kable", 
   add.median = TRUE,
   fixed.digits = TRUE, 
   simple.kable = FALSE,
   title="",
   numformat = NA) %>%
  kableExtra::kable_styling(latex_options = c("scale_down", "hold_position"),
                            full_width = FALSE, fixed_thead = T)
Variable N Mean Std. Dev. Min Pctl. 25 Pctl. 50 Pctl. 75 Max
i_1_p9 18995 5.761 1.363 1 5 6 7 7
i_2_p9 18995 5.513 1.326 1 5 6 7 7
i_3_p9 18995 5.135 1.922 1 4 6 7 7
i_4_p9 18995 5.213 1.554 1 4 5 6 7
i_5_p9 18995 5.926 1.284 1 5 6 7 7
i_6_p9 18995 4.855 2.057 1 4 5 7 7

¿Cuántos cluster?

Click para ver el código
# Importante: estandarizar
data_zscore <- scale(data)

fviz_nbclust(data_zscore, 
  kmeans, 
  k.max = 6,
  method = "wss") +  
  theme_minimal()

Aplicamos nuestra estrategia de clustering

Click para ver el código
set.seed(123)
model.km <- kmeans(data_zscore, 
                   iter.max=5,
                   centers=3, 
                   nstart=100, 
                   trace = FALSE)

summary(model.km)
             Length Class  Mode   
cluster      18995  -none- numeric
centers         18  -none- numeric
totss            1  -none- numeric
withinss         3  -none- numeric
tot.withinss     1  -none- numeric
betweenss        1  -none- numeric
size             3  -none- numeric
iter             1  -none- numeric
ifault           1  -none- numeric
Click para ver el código
# Veamos los cluster
data$cluster_kmedias <- model.km$cluster
head(data)
  i_1_p9 i_2_p9 i_3_p9 i_4_p9 i_5_p9 i_6_p9 cluster_kmedias
1      7      6      7      6      7      7               3
2      4      4      3      6      6      6               2
3      4      4      5      4      5      5               2
4      6      4      5      6      4      4               2
5      1      1      1      3      2      6               2
6      4      5      6      6      5      4               2

Exploremos las categorías obtenidas:

Click para ver el código
# Visualicemos los cluster 
fviz_cluster(model.km, data = data, geom = "point") +
  labs(title=NULL) +
  theme_light()

Click para ver el código
# Exploremos los cluster
aggregate(data[1:6], by=list(model.km$cluster), mean)
  Group.1   i_1_p9   i_2_p9   i_3_p9   i_4_p9   i_5_p9   i_6_p9
1       1 6.161019 5.594303 2.696252 4.794903 6.128936 2.057871
2       2 4.331965 4.116706 4.154528 3.824941 4.523882 3.974065
3       3 6.239917 6.073779 6.284148 5.921961 6.452280 6.070788

Ejercicio propuesto

Piense en una pregunta de investigación que pueda ser respondida con datos con de la encuesta ENSSEX y proponga de forma intuitiva un modelo para este fenómeno. Implemente este modelo con las herramientas vistas en el workshop.